Molecular dynamics of folding of secondary structures in Go-type models of proteins 
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We consider six different secondary structures of proteins and construct two types of Go-type off- 
lattice models: with the steric constraints and without. The basic aminoacid-aminoacid potential 
is Lennard Jones for the native contacts and a soft repulsion for the non-native contacts. The 
interactions are chosen to make the target secondary structure be the native state of the system. We 
provide a thorough equilibrium and kinetic characterization of the sequences through the molecular 
dynamics simulations with the Langevin noise. Models with the steric constraints are found to be 
better folders and to be more stable, especially in the case of the /3-structures. Phononic spectra for 
vibrations around the native states have low frequency gaps that correlate with the thermodynamic 
stability. Folding of the secondary structures proceeds through a well defined sequence of events. For 
instance, a-helices fold from the ends first. The closer to the native state, the faster establishment 
of the contacts. Increasing the system size deteriorates the folding characteristics. We study 
the folding times as a function of viscous friction and find a regime of moderate friction with 
the linear dependence. We also consider folding when one end of a structure is pinned which 
imitates instantaneous conditions when a protein is being synthesized. We find that, under such 
circumstances, folding of helices is faster and of the /3-sequences slower. 



I. INTRODUCTION 

Understanding of the statistical mechanics aspects of 
protein folding has been recently advanced through stud- 
ies of coarse grained models in which aminoacids are rep- 
resented by single beads. In particular, many valuable 
insights have been gainacLby considering such models on 
a lattice (see, e.g. ref.ErEl). These toy models have al- 
lowed one to relate the folding j)rocess to the sequence 
dependent, energy landscapes,Qi3 to identify the foldiiig 
pathwaysjj and to study the issues of designabilityJa'El 
They have also been used to demonstrate existence of 
tree like kinetic connectivities in the folding funnelt3 that 
can be rfepresented with the use of the disconnectivity 
graphs .Ell 

More realistic coarse grained models, however, require 
an off-lattice setting. Recently, there have been a num- 
ber of off-lattice studies, which focus on the kinetics of 
foldingJl3"ll3 as well as on the sequence Hpsign and deter- 
mination of the interaction potentials .EZlEj It has been 
also suggested that the geometry of the native structure 
of proteins itself, without a detailed information regard- 
ing the amino acid— sequences, plays a decisive role in 
the folding process.t^ The effective interactions between 
the beads in the coarse-grained models are difficult to 
derive from microscopic considerations and instead they 
may be chosen to reflect statistical properties of protein 
structures as collected in the protein data bank.E3tJ 

The advantage of the lattice models is that they allow 
for an enumeration of conformations, at least for short 
chains, and thus for identification of the native state and 
determination of equilibrium properties of the system. 
However, the dynamics of these models are not related 
to any Newton's equations (in the classical limit) since 
they have to be defined in terms of the discrete Monte 
Carlo steps made within a declared set of allowed moves. 



Molecular dynamics (MD) simulations are a natural 
tool to study models set in a continuum space indepen- 
dent of whether they are coarse-grained or fully atomic. 
A decade ago, MD simulations of the microscopic repre- 
sentation of a polypeptide chain could explore time scales 
which were around nanoseconds and thus orders of mag- 
nitude too short to study the full duration of a typical 
foldingHa Curreatly, microsecond time scales have be- 
come accessibleCJ (in a special purpose computer) but 
these feats allow one to monitor individual trajectories 
in very restricted regions of the phase space. A detailed 
characterization is still restricted to up to 10ns long time 
scaleso Thus such chemically realistic models cannot yet 
provide a sufficiently thorough equilibrium and kinetic 
characterizations of the system that are expected when 
setting up a model to be studied within the framework 
of statistical mechanics. 

Thus there is a need to study simplified coarse-grained 
continuum space models to understand the generic fea- 
tures of folding. These models must involve idealized po- 
tentials, for instance of the Lennard- Jones kind. These 
interactions may either be constrained to correspond to a 
target native state or they are not.E3 The targeting may 
be facilitated by augmenijing the model by an introduc- 
tion of steric constraints, tilEj which takes into account 
the properties of the peptide bonds, but it can also be 
accomplished without such constraints. The dynamics 
of the simplified Lennard-Jones models are usually stud- 
ied by the methods of MD. A novel and efficient variant 
of the MD tedinique has been recently proposed by He 
and Scheraga,L3 in which one focuses the evolution on the 
torsional degrees of freedom. Sorrie .Monte Carlo studies 
for those models are also available ll3'E3tj The applicabil- 
ity of the Monte Carlo methods to dynamical properties 
(as opposed to equillibrium) remains, however, untested. 

In this paper, we report on molecular dynamics stud- 



ies of possibly the simplest models with the native states 
defined by target copformations - the ofF-lattice versions 
of the Go models.cZl The interesting property of the Go 
model is that it essentially avoids the issue of the correct 
specification of the aminoacid-aminoacid interaction and 
yet it corresponds to a realistic conformation by defin- 
ing the couplings in terms of the target conformation. 
Specifically, we consider the Lennard-Jones interactions 
between the beads and make them attractive for native 
contacts (two non-contiguous beads form a contact if 
they do not exceed a certain cutoff distance) and repul- 
sive for non-native contacts. Such models are protein- 
like also in the sense that they minimize the structural 
frustration. The object of our studies here is to investi- 
gate how viable are such Lennard- Jones-Go models in the 
context of the kinetics of protein folding. Our motivation 
for performing these studies follows from the expectation 
that such models may play a role similar to that of the 
Ising spin models in representing properties of the more 
complicated real life magnetic systems. 

Recently, there have been several related studies of 
the coarse-grained off-lattice Go models which, however, 
asked different questions than in this paper and were 
not based on the Lennard-Jones potentials. The study 
by Zhou and Karplusj23 employed a square well poten- 
tial (also used before to analyze homopolymersEH) which 
leads to a simplified discrete MD treatment. The au- 
thors have studied possible scenarios (with or without 
long-lasting intermediates) of the folding kinetics in a 
three-helix-bundle-like protein model as a function of 
the strength of the non-native contacts relative to the 
strength of the native ones. This kind of discrete MD 
have been also used by Dokholyan, Buldyrev, Stanley 
and Shakhnovichll3 to identify a folding nucleus in a Go- 
type heteropolymer. (For a related-full atom study of this 
issue for the CI2 protein see ref.c3). Another study, by 
Hardin, Luthey-Schulten, and WolynesEiJ used a detailed 
backbone representation and implemented an associative 
memory Hamiltionian in which the contact potentials are 
set to be the Gaussian functions. This study was aimed 
at understanding the kinetics of the secondary structures 
formation from the perspective of the energy landscape 
picture. Our partiality to the Lennard-Jones potentials 
is of a twofold nature. First, these potentials are well 
established in simulations of liquids and -essentially - 
continuum time stable MD codes are available. Second, 
their overall distance dependence seems qualitatively cor- 
rect on a fundamental level. 

Real life conformations of proteins in the native state 
consist of interconnected secondacy structures such as 
a-helices, /?-hairpins (see, e.g. ref.L3), and /3-sheets. Un- 
derstanding of the kinetics of protein folding should be 
first accomplished at the level of the secondary struc- 
tures. This is, in fact, the task of this paper and we 
narrow it to two classes of the Go models: with and 
without the steric constraints and for each of these we 
consider six possible secondary structures with different 
numbers of monomers. The kinetics of folding are stud- 



ied by the standard techniques of MD with a Langevin 
noise which controls the temperature of the systems and, 
at the same time, mimics the interactions of the protein 
fragments with the molecules of water. 

Our method of constructing the off-lattice Go mod- 
els is outlined in Section IL In Section III, we determine 
the sizes of the native bgisips of the models through the 
shape distortion methodcS instead of making assump- 
tions about them, as it is done usually. In Section IV, 
we determine the phononic spectra of the models and 
discuss their relationship to the thermodynamic stabil- 
ity of the native structures. In Section V, we deter- 
mine the folding times, the thermodynamic stability and 
other characteristic parameters for the models without 
the steric constraints and find that the /3-structures have 
very low thermodynamic stability in this model. The 
folding times vary with the native conformation and the 
chain length. In particular, the a-helix folds faster than 
the /3-hairpin of the same length. The latter also folds 
much slower than a similarly sized /3-sheet with three 
strands. Experiments on folding indicate that the fold- 
ing time for the hajicpin is about 30 times slower that 
that for the a-helix.o Our results do not yield the same 
rate but, at least, they show that the /3-hairpin is the 
slower folder. A similar result has been obtained by 
Hardin, Luthey-Schulten, and Wolynes.EiJ Our studies of 
the three sizes of the helices confirm the general observa- 
tion that folding ppaperties deteriorate with the growth 
of the system size.E3 In Section VI, we repeat the MD 
studies for the models with the steric constraints and find 
that these constraints substantially improve the thermo- 
dynamic stability, especially for the /3-structures. At the 
same time, however, they raise the temperature of the 
onset of the glassy effects so the net result is that the a- 
helices remain good folders and the /3-structures remain 
bad folders but their foldability is improved. Notice, how- 
ever, that our models allow for good or adequate folding, 
depending on the system, without any additional dipo]a= 
dipole interactions as introduced by He and Scheraga.o 

In Section VII, we focus on the mechanisms of folding 
and discuss sequencing of events that takes place in fold- 
ing of the secondary structures. In particular, we investi- 
gate characteristic time separations between the folding 
steps. In general, the time separations may depend on 
details of a model but the sequencing is expected to be 
model independent, i.e. it should proceed as predicted by 
our Go-type models. In Section VIII, we study the depen- 
dence of our results on the strength of viscous damping 
in the Langevin noise and find that it affects the folding 
time but it does not affect any of the other characteristic 
parameters except when the friction is set unrealistically 
low. Finally, in Section IX, we consider the processes of 
folding occurring when one end of the structure is fixed 
in space which imitates conditions found when a protein 
is being synthesized. We find that fixing one end of 
a helix accelerates folding and improves overall folding 
characteristics. The folding is found to originate at the 
undamped end. On the other hand, folding of the (3 



structures becomes worse when one end is clamped. 

This paper sets the stage for studies of fuU proteins and 
kinetic interplays between various secondary structures. 
Results of such studies for the Go models of proteins will 
be presented in a separate paper. 



II. MODELS 

We represent a polypeptide in a simplified manner: by 
a chain of connected beads. The beads' positions corre- 
spond to the locations of the Ca atoms. When the system 
is in its native state, the beads have no kinetic energy and 
the locations of the Ca atoms can be obtained from the 
Protein Data Bank. However, the secondary structures 
that we study are idealized and are not targeted to a 
specific real structure. Instead, they are just meant to 
be constructed in a way which is very close to typical 
a-helices and /3-sheets found in the native states of real 
proteins. These secondary structures are stabilized pri- 
marily by the hydrogen bonds. The structures that we 
study are shown in Figure 1. These are: three a-helices, 
denoted as HIO, H16 and H24, two /3-hairpins, denoted 
as BIO and B16, and a /3-sheet, B15. The numbers in the 
labels indicate the numbers of beads in the structures. 

All of the bonds between two consecutive beads along 
the chain of the helices, i.e. the peptide bonds, are as- 
sumed to have the same length of do = 3.8A which is a 
typical real life value. As one proceeds along the axis of 
the helix (the z-coordinate) from one bead to another, 
the bead's azimuthal angle is rotated by 100° and the 
azimuthal length is displaced by 1.5 A which again corre- 
sponds to a typical geometry. 

Each of the /3-hairpins, BIO and B16, has two anti- 
parallel strands which are connected by a turn. In the 
B15 /3-sheet there are three strands and two turns. In 
the /3-structures, the strands are not straight lines but 
have a zigzagged geometry as shown in Figure 1. The 
bond lengths between two connected beads are equal to 
3.8A but the displacement along the strand direction (the 
z axis in Figure 1) is 3.5A. The distance between two 
opposite beads in two bonded strands, is set to 5A, which 
is roughly equal to the hydrogen bond's length. The turn 
region is constructed so that the bond length and the 
zigzag pattern match. 

The potentials of interactions between pairs of the 
beads are constructed in a way that ensures that the 
target structure coincides with the ground state of the 
system, i.e. with the native state. We pick the pair po- 
tentials to be of the Lennard-Joaes type and select the 
parameters in a Go-type fashioriED so that significant at- 
traction is associated with the native contacts and the 
non-native contacts are purely repulsive. In our model, 
we assume that a native contact is present if the distance 
between the two monomers in the designed structure is 
shorter than 7.5A. 
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FIG. 1. Stereographic projections of the target structures 
studied in this paper: three a-helices HIO, H16, H24, two 
/3-hairpins BIO, B16, and a /3-sheet B15. The xz and yz 
planar projections for HIO and BIO are shown on the right 
hand side. 

These Go-type couplings already stabilize the struc- 
tures under studies but that stabilization will be found 
here not to be sufhciently adequate, especially in the case 
of the hairpins. We thus consider two classes of models: 
with and without additional steric constraints. The steric 
constraints add extra stability and they take into account 
the directional character of the peptide bonds in a more 
realistic manner. 



A. Go-type model with no steric constraints 

This model is similar iit-spirit to that introduced by 
lori, Marinafi, and Parisit3 and to that studied by Li 
and Cieplak.LJ The conceptual difference between these 
two papers is that the former is not constructed in refer- 
ence to any predetermined target structure. For a con- 
formation defined by the set of position vectors {r^}, 
i = 1,2...N, the potential energy is assumed to take 
the following form 
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The first term represents rigidity of the backbone poten- 
tial, the second term corresponds to interactions in the 
native contacts and the third term to those in the non- 
native contacts. N denotes the number of residues. 

The backbone potential takes the form of the sum over 
harmonict2l and anharmonicES interactions 



N-1 

= - do? + fc2(rM+i - do)% (2) 

i=l 

where r'i^i+i — jr^ ~rj_[_i | is the distance between two con- 
secutive beads; do = 3.8A, fci = e and /c2 — lOOe, where 
e is the Lennard- Jones energy parameter corresponding 
to a native contact. 

The interaction between residues that form a native 
contact in the target structure is taken to be of the 
Lennard-Jones form: 
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where the sum is over ah pairs of residues i and j (but 
those which are immediate neighbors along the chain) 
which form the native contacts in the given target struc- 
ture. Tij = jr^ — r_, | is the monomer to monomer distance. 
The parameters cry are chosen in a way that each contact 
in the native structure is stabihzed at the minimum of 
the potential Essentially, cry 
the corresponding native contact's length 

Residues not forming the native contacts interact via 
a repulsive soft core potential. Our potential for non- 
native contacts, given below, differs from the model of 
lori et al. in that it falls to after some cut-off distance, 
dcut- The purpose of the cutoff is to minimize structural 
frustration and thus to improve foldability. 
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Here, cro = 2~^/^ • dcut- For distances shorter than dcut 
the potential is purely repulsive. For the ct-helices we 
chose dcut = {dij) ~ 5. 5 A which seems to be a natu- 
ral choice. For the /3-structures, however, such a choice 
of dcut leads to an instability of the native conformation 
due to a substantial energy degeneracy of the nearby con- 
formations. A sufficient extending of the cutoff removes 
the degeneracy and stabilizes the native state. We have 
found that an adequate choice for the /^-structures is to 
take dcut = 7.5A - the distance that defines what is a 
contact. 

Notice that in the target structures, all native con- 
tacts are optimized while the non-native contact poten- 
tials give no contribution to the energy. Thus the tar- 
get structure is indeed the native state of the system. 
The non-native interactions contribute only when mov- 
ing away from the target. 



B. Go- type model with the steric constraints 



In this case, the potential energy of the system is given 

by 

Ep({rj) = v^"" + y^^^ + y^o^ + + y^^, (6) 

where the first three terms are identical to those in equa- 
tion (|l|), whereas the last two terms correspond to the 
bond angle and dihedral angle potentials respectively. 
The bond angle, 9i, is defined as the angle between two 
successive vectors and v^+i, where = r^+i — r^. The 
dihedral angle, 0^, is the angle between two-vector prod- 
ucts Vi_i X Vi and x v^+i. Following ref.EJ, we use the 
following potentials for the bond and the dihedral angles 
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where ke = 20e/{radf, A = Oe and B ^ 0.2e^ur angle 
dependent potentials differ from those in ref.EJ, since in 
our case we take ^oi to be, in general, site-dependent and 
equal to the bond angles in the native targets. For the 
helices, however, ^oi Stve uniform and equal to 1.56157 
rad (89.4714 deg). For the /3-structures, 6'i's in the turn 
region are different from those in the strands. 

Introduction of the steric constraints to the model de- 
scribed by eq. (1) shifts the native state away from the 
target conformation because the target need not corre- 
spond to a minimum of the dihedral potentials. However, 
we have found that for our choice of the parameters A 
and B the true native states differ only little from the 
targets. The true native states are found by a multi- 
ple zero temperature quench procedure from low energy 
conformations generated by MD trajectories that start 
in the target conformation. The conformation distances 
from the native states to the targets (to be defined in 
Sec. Ill) never exceeded 0.05A. 

It should be pointed out that none of the models 
studied here can distinguish between the right- and left- 
handed helices. The ct-helices found in nature, however, 
are only right-handed. The model in which only right- 
handed helices are favored should include terms related 
to the chirality of the chain conformation. The chirality 
can be given by signlivi-i x • v^+i], where a nega- 
tive sign means a left-handed conformatioH-and a positive 
sign means a right-handed conformation. Ell 

In the following Sections, the Go- type sequences with- 
out the steric constraints will be denoted by the symbols 
associated with the targets: HIO, H16, H24, BIO, B16, 
B15. On the other hand, the sequences constructed with 
the additional steric constraints will be labeled as HIO, 
H16, H24, BIO, B16, B15. 



C. Dynamical equations and the thermostat 

The motion of the model secondary structures of pro- 
teins can be described by the Langevin equation 
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where r is a generahzed coordinate of a bead, m is the 
monomer's mass, Fc = —VrEp is the conformation force, 
7 is a friction coefficient and F is the random force which 
is introduced to balance the energy dissipation caused 
by friction. Both the friction and the random force rep- 
resent the effects of the solvent and they control the 
temperatureHZI F is assumed to be drawn from the Gaus- 
sian distribution with the standard variance related to 
temperature by 



(F(O)r(i)) = 2^kBT6{t), 



(10) 



where fc^ is the Boltzmann constant, T is temperature, 
t is time and 6(t) is the Dirac delta function. 

The Langevin equations are integrated using the fifth 
order predictor-corrector schemeH The friction and ran- 
dom force terms are included in the form of a noise per- 
turbing the Newtonian motion at each integration step. 
In the case of the model with the steric constraints, 
the forces associated with the angle-dependent potentials 
are calculated through a numerical determination of the 
derivatives of the potential. 

In the following, the temperature will be measured in 
the reduced units of e/ks- The integration time step 
is taken to be At — O.OOSt, where r is a characteristic 
time unit. At low values of friction, t coincides with the 
period of oscillations, t„i near the Lennard- Jones min- 
imum and is equal to y^ma^ /e, where a is a Van der 
Waals radius of the amino acid residues. The value of 
a is chosen to be equal to 5 A, andJ-his value is roughly 
equal to (cTy ). As estimated in ref.Eil, the typical values 
of TO is 3 X 10~^^(7 and e is of the order of Ikcal/mol, 
hence r is roughly equal to 3ps. The period Tm, how- 
ever, depends on the friction coefficient, 7, as discussed 
in ref.EJIlj. Most of our simulations are performed with 
7 = 10tot~^ which is 5 times higher than a standard 
choice in molecular dynamics studies of liquids. Higher 
values of 7 may be actually more realistic. For a single 
aminoacid in water, an. effective 7 has been argued to 
be even of order SOto/t.EJO However, the MD code be- 
comes unstable for too high 7's. The dependence of the 
folding characteristics on 7 will be discussed later in this 
study. 



III. DETERMINATION OF THE NATIVE BASIN 

One problem that is encountered when studying off- 
lattice models is providing a definition of what small dis- 
tortion away from the native state can still be considered 
as belonging to the native state. In other words - what 



is the definition of the native basin (which should not 
to be confused with the folding funnel)? Answering this 
question is important when calculating the folding time 
(is the system already in the native basin?) and when 
calculating the equilibrium stability of the basin. A need 
for delineation of the native basin does not arise in lattice 
models due to the discretization of the possible shapes - 
there, it is just one conformation. 

The delineation of the native basin can be accom- 
plished by considering the conformational distance, S, 
away from the native state. S is given by 
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where and rfj^^ are the monomer to monomer dis- 
tances in the given structure and in the native state re- 
spectively. The distance 6 is measured in A. 
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FIG. 2. The average root mean square distance to the na- 
tive state as a function of time when the sequence is placed 
in its native conformation. The figure is for sequences H16 
and B16 and for the reduced temperatures as indicated. The 
results are averaged over 500 trajectories which differ in the 
strings of the random noise. 

We estimate the characteristic size of the native basin 
by using the shape distortion method proposed in ref.E3. 
The method is based on monitoring the conformational 
distance to the native state as a function of time and it 
considers trajectories that originate from the native state. 
At a sufficiently large time scale the conformational dis- 
tance saturates below some critical temperature, Tc- The 
saturation value of the distance at this temperature, Sc, 
is used for the estimated native basin's size and the re- 
sults are found to be consistent with the quench-obtained 
sampling of the local energy minima of the system. At 
the same time, Tc has been founded to be a measure of 



the folding temperature Tf . Tf itself is defined as a tem- 
perature at which the equilibrium probability of finding 
the system in its native basin crosses 1/2. 

An illustration of the shape distortion method for se- 
quences H16 and B16 is shown in Figure 2. For each 

1 /2 

of these sequences the dependence of (iJ^)^ on time is 
shown for two temperatures: a critical Tc, when the sat- 
uration is still observed, and a somewhat higher T, when 
there is no saturation which is interpreted as exiting the 
"trap" provided by the native basin. The estimated basin 
size for sequences H16 and B16 are 0.6A and O.ISA re- 
spectively. Notice that sequence H16 has much larger 
basin than sequence B16 and the corresponding Tc is also 
much higher for H16. 

Table I summarizes the values of Sc and Tc for all of 
the sequences studied. It also shows the number of native 
contacts in the native state, Nc, the conformational dis- 
tance from the native state to the closest local minimum, 
6i, and the values of Tf obtained in the next Sections. 
The closest local minimum to the native state is found by 
a multiple quenching from the random conformations in 
a low T MD trajectory. Notice that for the /3-structures, 
the models with the steric constraints yield a larger 6c 
compared to the models without the constraints. The 
basin sizes for the a-helices, on the other hand, are com- 
parable. In most cases, we have found that Sc is some- 
what smaller than Si. This is not so, however, in the 
case of H24 and B15. We interpret this as a necessity 
to delineate the native basin by more than one isotropic 
parameter Sc and yet we shall use it here as a simple way 
to provide estimates. 



SEQ 




SilA] 


Sc[A] 




Tf 


Tmin 


HIO 


21 


0.8248 


0.50(3) 


0.23(1) 


0.24(1) 


0.25(3) 


H16 


39 


0.6594 


0.60(3) 


0.24(1) 


0.24(1) 


0.30(3) 


H24 


63 


0.5430 


0.53(3) 


0.19(1) 


0.20(1) 


0.30(3) 


BIO 


18 


0.2588 


0.15(2) 


0.035(2) 


0.033(2) 


0.08(1) 


B16 


33 


0.2288 


0.15(2) 


0.030(2) 


0.027(2) 


0.07(1) 


B15 


33 


0.1382 


0.13(2) 


0.030(2) 


0.029(2) 


0.05(1) 


HIO 


21 


0.8133 


0.46(3) 


0.30(1) 


0.31(1) 


0.50(5) 


H16 


39 


0.6356 


0.55(3) 


0.30(1) 


0.29(1) 


0.45(3) 


H24 


63 


0.3896 


0.55(5) 


0.24(2) 


0.24(1) 


0.36(3) 


BIO 


18 


0.3178 


0.25(2) 


0.20(1) 


0.20(1) 


0.35(5) 


B16 


33 


0.2801 


0.25(2) 


0.20(1) 


0.21(1) 


0.35(2) 


B15 


33 


0.2045 


0.30(2) 


0.22(1) 


0.20(1) 


0.35(2) 



TABLE I. The number of native contacts in the native con- 
formations, Nc, the conformational distance from the native 
state to the closest local minimum, 5\, the native basin size, 
5c, the critical temperature of the shape distortion, Tc, the 
folding temperature, T/, and the temperature of the fastest 
folding, Tmin, for the sequences studied. The numbers in 
parenthesis indicate the error bars. 



IV. PHONONIC SPECTRA AND STABILITY 

We now consider elastic vibrations of the systems 
around their native states and ask if the phononic spectra 
relate to the folding temperature. 

Let Uia denote a small displacement of bead i in the 
direction a [a — x, y, z) from its native position. In 
the harmonic approximation, the motion oLthe system is 
governed by a set of 3iV coupled equationsLJ 
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where lo is an angular frequency, and k"^^ are the second 
derivatives of the potential energy taken at the native 
conformation: 
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We diagonalize these equations with the use of the Jacobi 
transformation methodcj and determine the phononic 
spectrum from the eigenvalues. The elastic constant ma- 
trix {kfj^} is real and symmetric. For the models without 
the steric constraints, when only the pair-wise interac- 
tions are present, the diagonal elements are equal to the 
negative of the sum over all of the off-diagonal elements 
in the same row (or column). In the case of the models 
with the steric constraints, however, the potentials ac- 
quire many-body components. The elastic constants k°j^ 
are calculated by freezing all beads but one, at a time, 
in their native positions and by measuring the resulting 
increase in the force when the unfrozen bead is displaced 
in a given direction. This is done numerically and then a 
displacement of lO^^A has been found to yield a sufhcient 
accuracy. The elastic constants, when measured along a 
line connecting the two beads, are of order of 2e/A^ for 
both the contact and peptide bond interactions. The 
steric constraints enhance the elastic constants by up to 
a factor of 2. Naturally, this will increase the stability. 

The phononic spectra for the models without and with 
the steric constraints are shown in Figures 3 and 4 respec- 
tively. Each spectrum contains 6 zero frequency modes 
which correspond to the translational and rotational de- 
grees of freedom of the system as a whole. The first ex- 
cited mode, of frequency wi, defines the frequency gap in 
the spectrum. We observe that the gap correlates with 
the thermodynamic stability: generally, the bigger the 
Tf, the bigger the gap. Specifically, we notice that, for a 
given kind of structure, uji decreases as the system size 
increases. Furthermore, the gaps for the a-helices are 
larger than those for the /3-structures for a given chain 
length. Finally, adding the steric constraints shifts the 
whole spectra towards higher frequencies. 

An energy associated with a mode depends on its am- 
plitude. If an amplitude of vibrations of a single bead is 
of order uq, then the energy of the first excited mode, Ei, 



is of order NitiujIuq. As an estimate of uq, we may take 
« 0.1a which is an analog of the Lindemann criterion for 
melting of the solids. Note that uji is expressed in units 
of 1/r, where t = ^JmaFJe and a = hA. Thus our esti- 
mates for E-yje are 1.32, 0.77 and 0.28 for HIO, H16 and 
H24 respectively. This indicates a decrease of stability 
with N but the correlation with the actual value of Ty is 
weak. 
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FIG. 3. The phononic spectra for the sequences modeled 
without the steric constraints. The values of wir given in rad 
are, top to bottom: 3.635, 2.195, 1.080, 1.789, 0.815, 1.339. 
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FIG. 4. Same as Fig. 3 but for the sequences modeled with 
the steric constraints. The values of uj\t are, top to bottom: 
4.385, 2.290, 1.105, 2.775, 1.695, 1.795. 



V. FOLDING PROPERTIES IN MODELS 
WITHOUT THE STERIC CONSTRAINTS 

A. Equilibrium properties 

Before we discuss the folding process, we first provide 
a thermodynamic characterization of the sequences stud- 
ied. We focus on three parameters: Vq, C, and x- These 
denote the probability of being in the native basin, the 
specific heat per bead, and the structural susceptibility 
respectively. The thermodynamic stability may be char- 
acterized by T/ at which T'o crosses \. The specific heat 
is defined by the energy fiuctuations: 



C = 



1 (i?') - ^Ef 



N 



2^2 



(14) 



where E is the total energy (kinetic and potential) of the 
system. The brackets denote the thermodynamic aver- 
age. The structural susceptibility is defined in terms of 
the structural overlap functioiJlil which is given by 



N-2 N 



Xs = l 



1=1 j=i+2 



NAT 



(15) 



where Q{x) is the Heavyside function, x is then a mea- 
sure of fluctuations in Xs- 



X^{xl{T))-{Xs{T)Y 



(16) 



The maximum in C, when plotted against T, corre- 
sponds to the collapse transition temperature Tq at which 
there is a transition from random coil to compact con- 
formation. The maximum in x niay be interpreted as a 
signature of the |£ejlding temperature Tf .oc3 It has been 
suggested in ref.E^ that a small difference between these 
two temperatures is indicative of good folding properties. 
As a practical criterion for good foldability one may ts 
the parameter ctt = (Te - Tf)/Te to be less than 0.4.lli 

We calculate the thermodynamic parameters by aver- 
aging over many long MD trajectories using the native 
state as the starting configuration to maker sure that the 
evolution take space in the right parts of the phase space. 
For each temperature, the times used for averaging in 
each trajectory are between 500 and 2000t depending on 
the chain length. The first lOOOr are not taken into ac- 
count when averaging. We used as many trajectories as 
needed to provide stable characterization. In practice, 
considering of order 50 trajectories have been found to 
be sufficient. 

The values of Tf are listed in Table I. They almost co- 
incide with the values Tc obtained by the shape distortion 
method. One can clearly see that the stability of the he- 
lices is substantially larger than that of the /3 structures. 
The poor stability of the /3-systems is a result of their 
flat non-compact native conformation. There are many 



competing local energy minima with very similar struc- 
tures in the immediate neighborhood of the native state. 
We have checked that an increase in the off-plane zigzag 
like departures can bring about a significant boost to the 
stability. However, one of the goals of this paper is to 
deal with realistic geometries. We shall see, in Sec. VI, 
that the extra stability that is needed will be provided by 
the additional steric constraints. This is also illustrated 
in Figure 5 for three selected sequences. 
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FIG. 5. The probability of being in the native basin as 
function of temperature for the sequences with 7V=16. H16 
and B16 are without the steric constraints. Adding the steric 
constraints improves the stability as seen in the case of B16 
on the right hand side. 
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FIG. 6. The specific heat, C, and the structural suscepti- 
bility, X) as functions of the reduced temperature for the se- 
quences indicated. The arrows show the values of the folding 
temperature as obtained through monitoring of the equilib- 
rium probability to reside in the native basin. 

Figure 6 summarizes the thermodynamic properties 



of the sequences modeled without the steric constraints. 
The peaks of C and x are observed for all of the he- 
lices and /3-systems. Notice, however, that in the case 
of the /3-structures the peaks are much less pronounced. 
Notice also that the values of Tf , as indicated by the ar- 
rows and obtained based on monitoring Vq{T), generally 
agree quite closely with the positions of the maxima in 
X- This consistency testifies about a reasonable and con- 
sistent determination of the size of the native basin. For 
all of the sequences studied here the parameters ax are 
smaller than 0.4 and thus are expected to be good folders 
whereas H24 shows the borderline behavior. 



B. Kinetic properties 

We now consider the kinetic properties of the sequences 
without the steric constrains. From a kinetic point of 
view, a sequence is considered to be a good folder if 
its folding temperature Tf is comparable or preferablely 
larger than Tmin ~ the temperature of the fastest fold- 
ing. If this happens then there must exist a temperature 
range below Tf at which the native state is still easily 
accessible kinetically. Thus the first task in this context 
is to determine Tmin- 

We have studied the dynamics of the sequences by 
extensive MD simulations. Figure 7 shows the median 
folding time as a function of temperature. At each tem- 
perature the median folding time is computed based on 
either 100 or 200 trajectories starting from random initial 
conformations. The initial conformations are generated 
from random sets of the bond angles and dihedral angles, 
where the bond angles are additionally restricted to have 
random value between 0° and 90° so that these confor- 
mations are not compact. On the other hand, the cases 
with the distance between two non-bonded beads of less 
than O.&do have been also excluded in order not to violate 
the self- avoidance condition. Conformations generated 
by this method typically look like extended random coils. 
Figure 7 shows the expected U-shape dependences of the 
folding time on temperature. Notice that the /3-sequences 
fold the best at a much lower range of temperatures than 
what is observed for the a-helices. The temperature of 
the fastest folding Tmin for the /?-sequences is between 
0.05 and 0.08 while for the a-helices it is about 0.25 or 
0.3. 

Since for all the sequences Tf is less than Tmm, the se- 
quences can be considered to be either good or bad folders 
depending on how significantly do the folding conditions 
deteriorate on moving away from Tmin to Tf. Thus it 
is of interest to compare the folding time at Tmim ti, to 
the folding time at Tf, t2- The a-helices are good folders 
because for them ^2/^1 < 2. i.e., the two times are almost 
the same. For the the /3-sequences, on the other hand, 
these ratios are higher reflecting a much narrower tem- 
perature range in which the good folding conditions are 
available. For BIO and B15, for instance, the ratios t2/ii 



are about 7 and 5 respectively. Thus the /3-sequences 
without the steric constraints are bad folders. Note that 
our kinetically derived results on the quality of folding 
do not fully agree with those obtained by studying the 
thermodynamically derived parameter gt- The reason 
why the thermodynamic criterion appears to be faulty 
may have to do with the fact that the /3-structures stud- 
ied here have a low degree of compactness. Note also the 
poor pronunciation of the peaks in the specific heat for 
the /3-structures. 
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FIG. 7. The median folding time versus temperature for 
the sequences modeled without the steric constraints. The 
asterisks indicate the corresponding values of Tf. 

Figure 7 shows that the /3-sequences fold more slowly 
than the a-helices of the same chain length. For instance, 
the minimum folding time for /3-hairpin B16 is about 5 
times larger than for the a-helix H16. Notice also that 
the /J-sheet B15 folds considerablely faster than the hair- 
pin B16, though their lengths are comparable. All of 
these observations are consistent with the large role of 
the local contacts (as measured along the chain) in es- 
tablishing good foldabilityO 

Figure 7 indicates that the folding time strongly de- 
pends on the chaiH-|length. Studies of Jattice models, 
first by Gutin et alsB and recently by us,E3 have pointed 
out that the folding time grows as a power law with the 
system size t ^ N^. Our calculationsE3 have shown that 
for the Go sequences A is about 6 for two dimensional 
lattice sequences and about 3 for three dimensional se- 
quences. Here, we have checked that the folding times at 
Tmin of the three sequences HIO, H16, H24 are consis- 
tent by the power law with A — 3.0 ±0.1. The agreement 
with the exponent obtained for the maximally compact 
Go lattice models is probably coincidental. One of the 
important conclusions suggested in refo is that there 



may exist a limit to functionality of proteins that can be 
demonstrated through studies of scaling of Tf and Tmin- 
In the Go lattice models, as TV grows, Tmin grows indefi- 
nitely whereas Tf saturates with TV. This, together with 
the increase in the folding times, indicates a deterioration 
of the folding properties with N. The data obtained for 
the few sizes of helices studied here (see Table I) are only 
partially consistent with this scenario for the behaviors 
of the characteristic temperatures. Note, however, that 
the longer a helix, or a /J-structure, is the less compact 
is its native conformation and this in itself may yield a 
different behavior of Tf and Tmin- 



VI. FOLDING PROPERTIES IN MODELS WITH 
THE STERIC CONSTRAINTS 

A. Equilibrium properties 

We now repeat our methodology in reference to the 
sequences which are modeled with the implementation 
of the steric constraints. The thermodynamic properties 
for the systems of this type are summarized in Figure 8. 
The values of the folding temperatures Tf are also listed 
in Table I. Notice that the models with the steric con- 
straints are generally characterized by larger values of Tf 
than the models without such constraints. The changes 
are especially substantial for the /3-structures since the 
resulting values of Tf become now comparable to those 
of the helices. 
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FIG. 8. Same as Figure 6 but for the sequences modeled 
with the steric constraints. Notice the shift in the tempera- 
ture scale. 



B. Kinetic properties 



A. Sequencing of events 



Figure 9 shows the dependence of the folding time on 
temperature. The points shown are based on either 100 
or 200 trajectories starting from random initial confor- 
mations. For all of the structures studied here, T^in 
also becomes larger than that in the models without the 
steric constraints. The time scales of folding, however, 
are affected to a relatively much smaller extent. For the 
/3-sequences, the folding time at T^m practically does not 
change. For the a-helices, on the other hand, some slow- 
ing down is observed. When comparing the a-helices to 
the /3-structures one concludes that the helices continue 
to be faster folders but, for instance, H16 now folds only 
about two times faster than B16 at T„i„, as opposed to 
the factor of five observed without the constraints. 

Like in the model without steric constraints, for all the 
cases we still have Tf < Tmin- Comparing the folding 
times at Tf and at Tmin yields that the helices continue 
to be good folders {t2/ti < 2). The /3-sequences are 
still bad folders but now the ratios t2/ti are somewhat 
smaller than that in the model without steric constraints. 
For instance, for BIO, ^2/^1 is about 3. 

We conclude that adding the steric constraints to the 
Go-type interactions between the aminoacids improves 
the stability of the sequences, while it does not change 
much in the folding characteristics. The folding charac- 
teristics seem to be embedded in the very geometry of 
the native structure. 
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FIG. 9. Same as Figure 7 but for the sequences modeled 
with the steric constraints. 



In this Section, we focus on identifying stages in the 
folding process of the model secondary structures. Specif- 
ically, we discuss orders and time scales in which specific 
contacts are established. We narrow our discussion now 
to the three following structures: the a-helix H16, the 
/?-hairpin B16, and the /?-shcct B15. 

By starting from random unfolded conformations wc 
compute an average time, at which a given contact is 
established for the first time. Two monomers, i-th and 
j-th, are assumed to be in a contact if their distance is 
smaller than l.bcfij. The contact range definition of l.bA, 
used in the potential design procedure, is not appropriate 
to be used here, since each contact has a potential with 
its own variable Uij . Since the number of all native con- 
tacts is quite large we select several contacts, as indicated 
in Figure 10, for monitoring. The choice of the contacts 
is motivated by their dominant role in stabilization pro- 
vided by the hydrogen bonds. The labels indicate order 
in which these selected contacts first appear on average. 
For the a-helix, the folding is seen to start, statistically, 
at the ends (we have symmetrized the data with respect 
to the two ends). On the other hand both /3-structures 
first fold near the turns, i.e. where the local contacts are 
present. Note that the contacts indicated in the helix all 
have the same locality index and yet the end contacts are 
found to be privileged. Adding the steric constraints does 
not change the average sequencing of the kinetic events. 

Figures 10 and 11 show the average times to estab- 
lish the contacts at Tmin for the unconstrained and con- 
strained Go-type models respectively. When comparing 
the times obtained for the a-helix to the /3-structures, 
note the change in the time scale. The full folding takes 
place in about 10 to 20 times longer than the time needed 
to establish all contacts. The general pattern observed 
here is that of steps which are initially more or less evenly 
spaced followed by a rapid acceleration towards the end. 



VII. THE MECHANISM OF FOLDING OF 
SECONDARY STRUCTURES 




H16 


T=0.3 

12 3 4 567 









15 <t>/T 


B16 


T=0.05 




1 


2 3 4 5 6 7 





B15 


T=0.06 
1 2 3 4 


150 <t>/T 



150 <t>/T 



FIG. 10. Top: the contacts considered in the sequencing 
studies together with the number which indicates the order of 
the appearance. The structures are, from the left to the right, 
H16, B16 and H15. Bottom: the corresponding average time 
for a contact to appear for the first time during folding. The 
folding takes place at Tmin- The results for H16 are averaged 
over 1000 trajectories that start from random conformations. 
For B16 and B15, 500 trajectories are considered. The error 
bars are indicated by the thicker horizontal marks. 
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FIG. 11. Same as the bottom of Figure 10 but for the ster- 
ically constrained sequences. 



Figure 12 shows the average times to establish the con- 
tacts at temperatures which are away from Tmin - above 
or below. We focus on H16 and B16. The results show 
that for sequence B16 the sequencing of contacts does not 
change with temperature. For H16, on the other hand, 
the sequencing starts the same as at Tmim i-e. at the 
end points, but the canonical order becomes disrupted 
at the later stages. For instance, establishing contact 5 



before 7 is as likely as establishing 7 before 5. Notice 

that, above T„ii„, establishing the basic contacts in the 
two structures takes about as long as at Tmin even though 
the folding times got longer. On the other hand, below 
Tm,ini this process becomes more extended in time. 
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FIG. 12. The contact sequencing for sequences H16 and 
B16 at temperatures that are above or below Tmin- 

B. Examples of the trajectories 

Notice that the first appearance of the last contact that 
is monitored docs not coincide yet with the full folding 
because a substantial time is still needed to lock precisely 
into the native basin. In fact, getting to the stage where 
the last contact starts giving contribution to the energy 
takes of order of only from 5% to 10% of the full folding 
time. This time roughly corresponds to the collapse of 
the chain. 

This is illustrated in Figures 13 and 14, for H16 and 
B16 respectively, where examples of single trajectories at 
Tmin are shown. The trajectories are characterized by 
their total energy, E, potential energy, E^, the radius of 
gyration, Rg, conformational distance away from the na- 
tive state and the number of all native contacts, iVg, (not 
only of those few that were discussed in the context of 
sequencing) present. All of these parameters depend on 
time essentially monotonically. The establishment of the 
short contacts takes place rapidly and then the system 
spends most of the folding time "floating" near the native 
basin until it really finds it. Thus what takes most of the 
time in folding is the process of search. This process is 
thus governed more by kinetics than by energetics. 
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FIG. 13. A typical folding trajectory at Tmin for sequence 
H16. In the top panel the left arrow indicates when the last 
contact appears for the first time, and the right arrow indi- 
cates when folding takes place, i.e. S < Sc. The horizontal 
dash line indicates the energy of the native state. 
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C. Dependence on the viscous friction 

We now examine how do the folding properties depend 
on the viscosity. So far, the results of all of the simula- 
tions shown here were obtained with 7 — 10m/ t. We 
have checked that for sequence H16, T„ii„ does not de- 
pend on the friction and is equal to 0.3 if 7 > 0.2m/T. 
At very low friction, i.e. when 7 < 0.2m/ t, we observe 
some decrease of Tmin. For instance, Tmin is about 0.05 
for 7 — 0.005m/ T. Since Tf cannot depend on 7, at this 
very low value of friction the sequence becomes an ex- 
cellent folder. Howeverj-such small values of 7 are very 
unrealistic for proteins .O'tj We have also checked that 
other thermodynamic properties do not depend on 7. 
Figure 15 shows the dependence of the median folding 
time on friction for a wide range of 7, at T = 0.3 and 
also at Tmin for 7 < 0.2m/r. Notice that there is a min- 
imum of the folding time at some intermediate friction. 
The folding time grows not only at higher friction, as 7 
increases, but also at very low friction, as 7 decreases 
towards 0. There is no folding at 7 = 0, when the total 
energy of the system is conserved - a coupling to a heat 
reservoir is essential for folding. The optimal value of 
7, at which folding is the fastest, is in a range of from 
0.2 to 0.5. This observation is in agreement with the 
results of Klimov and Thirumalai on the viscosity de- 
pendence of the folding rate.E3 Notice that Klimov and 
Thirumalai have determined foldability rates and not the 
folding times themselves. We also observe that for 7 be- 
tween 1 and 30, the folding time of sequence H16 grows 
almost linearly with 7, as shown in the inset of Figure 
15. By extrapolation we estimate that for 7 — 50m /t the 
folding time for H16 would be ISOOr (for H16 it would 
be about 3000t). Veitshans et alt3 have estimated that 
at this apparently more realistic value of friction, corre- 
sponding to the case of an single aminoacid in a water 
solution at room temperature, the time unit of r should 
be about 3ns. Hence our estimate of the folding time of 
a 16-monomer a-helix is from 4 to 10 /is. If we assume 
that for the /3-structures the folding time also scales lin- 
early with 7, then the time scale for the formation of 
a 16-monomer /3-hairpin would range from 24 to 30 us. 
For comparison, experimental results have shown that 
a /3-hairpin consisting of 16 aiainoacids can fold as fast 
as 6^3 at room temperature.Eil An a-helix of the same 
length folds 30 times faster. Notice that our estimates 
yield larger folding times than observed in experiments 
but the orders of magnitude agree. 



FIG. 14. Same as Figure 12 but for sequence B16. 
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FIG. 15. The median folding time at T — 0.3 (continuous 
line) as function of the friction coefflcient, 7, for sequence 
H16. For 7 < 0.2m/T, Tmin becomes smaller than 0.3 and 
the folding times at T = Tmin are shown for several values of 
7 (discontinuous line). For 7 > 0.2m/r, Tmin is equal to 0.3. 
The points are based on 200 trajectories. The inset shows the 
dependence at T = 0.3 in a non-logarithmic scale of 7, where 
for the higher values of 7, the points are fitted by a straight 
line. 



D. Folding with one end fixed 

Finally, we examine a special case, when one end of the 
chain is fixed during folding. This is relevant to the pro- 
cess of protein synthesis: one end can be thought of as 
being momentarily glued to the surface of a ribosome 
until the jjarhole protein is constructed by adding new 
segments.L^ In the synthesis process, the protein folds as 
it is being produced and its length extends. At each in- 
stant, however, one end of the protein can be considered 
pinned. 

Figure 16 shows that the helix H16 with the end bead 
fixed has a considerably lower value of Tmin, compared to 
no clamping, and it generally folds faster at low tempera- 
tures, pinned folds faster below Tmin and is characterized 
by a significantly lowered Tmin ■ We have also found that 
Tf is not affected by the pinning. Thus the helix be- 
comes a much better folder. On the other hand, as seen 
in Figure 16, the contact sequencing becomes somewhat 
disturbed. Folding starts at the undamped end. After 
the initial stage, many contacts get established almost 
simultaneously in an avalanche-like process. Then there 
is a large gap in time before the last contact, near the 
clamped end, comes into place. 

Analyzing the phononic spectrum of the pinned helix 
H16 yields loit « 1.964 for the first excited mode which 
is almost the same as for H16 without the pinning. 



For the /^-sequences the different scenario is different: 
the pinned /3-structures seem to fold slower than the un- 
pinned ones. For instance, fixing one end of sequence B16 
increase the folding time by about 60%, at T = 0.05. The 
pinned BIO also folds slower, but not significantly slower. 

Our results indicate that the process of protein synthe- 
sis itself may accelerate the folding process when a helix 
part is being produced. Short /3-structures may introduce 
some slowing down. 
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FIG. 16. The order of contact appearance together with 
the average times for the contacts to be established for the 
first time for sequence H16 with one end monomer fixed. The 
fixed monomer is shown enlarged. The results are averaged 
over 1000 trajectories at Tmin- The inset compares the tem- 
perature dependencies of the median folding time for the sit- 
uation in which the end monomer is (solid line) or is not 
(broken line) clamped. 



VIII. CONCLUSIONS 

In summary, we have studied the Go-type models of 
off-lattice secondary structures of proteins. We have pro- 
vided a systematic characterization of both equilibrium 
and kinetic properties of these models. Models with the 
steric constraints were endowed with better thermody- 
namic stability. The stability can be assessed from the 
phononic spectra. The folding times strongly depend on 
the system size and on the geometry of the native state. 
The a-helices, which are stabilized merely by the local 
contacts appear to have much better folding properties 
than the /3-structures. However, it should be noted that 
for both kinds of the structures none of the models is 
a very good folder - Tf is always smaller than Tmin, ~ 
although the choice of the native basin size Sc has been 
made consistently by the shape distortion technique. We 
have checked that if one uses a somewhat larger 5c then 
this results in an increase in the effective values of both 



Tf and T„iin and thus their ordering remains unchanged. 
The kinetic criterion for good foldabihty does not appear 
to be compatible with the thermodynamic criterion pro- 
vided by Khmov and ThirumalaCj in the case of the (3- 
sequences. We speculate that this may be related to the 
low compactness level of their native structures. Folding 
of the secondary structures at Tmin proceeds, on aver- 
age, through a well defined sequence of events. That se- 
quencing may become more complicated in proteins when 
several secondary structures compete when folding. We 
intend to explore this issue later. 
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